library(mvtnorm)
library(copula)
library(energy)
library(fGarch)
simulation.runs = 1000
n = 1000

x1 = rnorm(n)
p.x2x1=0.9
p.yx1=0.9
x2=sqrt((1-p.x2x1)^2)*rnorm(n)+x1*p.x2x1
y=sqrt((1-p.yx1)^2)*rnorm(n)+x1*p.yx1
y.reg=lm(y~x1)
e1=y.reg$residuals
x2.reg=lm(x2~x1)
e2=x2.reg$residuals
cov(cbind(y,x1,x2))
p.value = dcov.test(e1,e2,R=199)$p.value
plot(e1,e2)


###Multivariate T
df=5
cov.x2x1=.6
cov.yx1=.6
sigma=diag(3)
sigma[1,2]=cov.yx1
sigma[2,1]=cov.yx1
sigma[3,1]=cov.x2x1
sigma[1,3]=cov.x2x1
multivariate.t = rmvt(n,sigma,df)

pairs(multivariate.t)
cor(multivariate.t,method="pearson")

y=multivariate.t[,1]
x1=multivariate.t[,2]
x2=multivariate.t[,3]
hist(y)
hist(x1)
hist(x2)
y.reg=lm(y~x1)
e1=y.reg$residuals
x2.reg=lm(x2~x1)
e2=x2.reg$residuals
cov(cbind(y,x1,x2))
p.value = dcov.test(e1,e2,R=199)$p.value
plot(e1,e2)
cor(e1,e2)
qqplot(e1,e2)


#multi-variate normal marginals, t-copula
n=1000
p.x2x1=.7
p.yx1=.6
p.yx2=0
nu=1.5
param=c(p.yx1,p.yx2,p.x2x1)
t.copula=tCopula(param,dim=3,dispstr="un",df=3)
normal.copula=normalCopula(param,dim=3,dispstr="un")

marginals = rep("t",3)
#marginals = rep("std",3)
#marginals = rep("norm",3)

marginals.params = list(list( df = nu),list( df = nu),list( df = nu))

#marginals.params = list(list(mean = 0, sd =1, nu = nu),list(mean = 0, sd =1, nu = nu),list(mean = 0, sd =1, nu = nu))
#marginals.params = list(list(mean = 0, sd =1),list(mean = 0, sd =1),list(mean = 0, sd =1))

#mvariate.with.copula = mvdc(t.copula,marginals,marginals.params)
mvariate.with.copula = mvdc(normal.copula,marginals,marginals.params)

random.multivariate.with.copula = rmvdc(mvariate.with.copula,n)
y=random.multivariate.with.copula[,1]
x1=random.multivariate.with.copula[,2]
x2=random.multivariate.with.copula[,3]
pairs(random.multivariate.with.copula)
y.reg=lm(y~x1)
e1=y.reg$residuals
x2.reg=lm(x2~x1)
e2=x2.reg$residuals
cov(cbind(y,x1,x2))
p.value = dcov.test(e1,e2,R=199)$p.value

plot(e1,e2)
cor(e1,e2,method="pearson")
qqplot(e1,e2)

######
n=100
mean=0
sd=1
beta.1=.8
beta.2=.7
nu=2.1
rate=1
#x1=rstd(n,mean=mean,sd=sd,nu=nu)
#x1=rt(n,df=nu)
x1=rexp(n,rate=rate)
#x1=rnorm(n)

e1.model=rnorm(n,mean=mean,sd=sd)
e2.model=rnorm(n,mean=mean,sd=sd)

#e1.model=rt(n,df=nu)
#e2.model=rt(n,df=nu)



y=beta.1*x1+e1.model
x2=beta.2*x1+e2.model

y.reg=lm(y~x1)
e1=y.reg$residuals
x2.reg=lm(x2~x1)
e2=x2.reg$residuals
cov(cbind(y,x1,x2))
cor(cbind(y,x1,x2))
pairs(cbind(y,x1,x2))

p.value1 = dcov.test(e1,e2,R=199)$p.value
p.value1
p.value2 = dcov.test(y,x2,R=199)$p.value
plot(e1,e2)
#cor.test(y,x2,method="pearson")
cor.test(e1,e2,method="pearson")
qqplot(e1,e2)
